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ABSTRACT 

Exoplanets will remain spatially unresolved for the foreseeable future. It is in principle 
possible to determine a planet's brightness markings by analyzing its disk-integrated 
brightness variations, in either thermal or reflected light. In order to identify the de- 
generacies inherent in these inverse problems, we develop analytic solutions to the 
associated forward problems. Our approach is to compute the "harmonic lightcurves" 
resulting from spherical harmonic maps of intensity or albedo. The convolutions of- 
ten contain a nuUspace: a class of non-zero maps that have no lightcurve signature. 
We derive harmonic thermal lightcurves for both equatorial and inclined observers. 
The nuUspace for these two situations is significantly different, with odd modes be- 
ing present in the latter case, but not the former. We therefore suggest that the 
Fourier spectrum of a thermal lightcurve is sufficient to determine the orbital incli- 
nation of non-transiting short-period planets, the rotational inclination of stars and 
brown dwarfs, and the obliquity of directly imaged planets. In the best-case scenario of 
a nearly edge-on rotator, factor-of-two measurements of the amplitudes of odd modes 
in the thermal lightcurve provide an inclination estimate good to a few degrees. In 
general, however, inclination estimates will remain qualitative until detailed hydrody- 
namic simulations and/or occultation maps can be used as a calibrator. We further 
derive harmonic reflected lightcurves for tidally-locked planets; these are higher-order 
versions of the well-known Lambert phase curve. We show that a non-uniform diffusely- 
refiecting planet with a precisely Lambertian phase curve may have planetary and 
Bond albedos significantly different from that inferred if the planet is assumed to be 
uniform. Lastly, we provide low-order analytic expressions for harmonic lightcurves 
that can be used for fitting observed photometry; as a general rule, edge-on solutions 
cannot simply be scaled by sini to mimic inclined lightcurves. 



1 INTRODUCTION 
1.1 Motivation 



Extrasolar planets are sufficiently small and distant that they will remain spatially unresolved for the foreseeable future. It is 
nonetheless possible to infer spatial inhomogeneities on these bodies through a) occultations, as when a planet passes behind 
its host star (Majeau et al. 2012; de Wit et al. 2012), or b) orbital and rotational motion (e.g., Knutson et al. 2007; Cowan 
et al. 2009). At thermal wavelengths, such maps probe the magnetic fields of stars, large-scale atmospheric circulation of 
sub-stellar objects, as well as thermal inertia and albedo inhomogeneities of airless bodies. Exo-cartography with reflected 
light, on the other hand, constrains the albedo markings of planets. 

In the current study we consider photometric variability due to rotational and orbital motion, i.e. method b). We seek an 
analytic solution to the convolution relating intrinsic spatial inhomogeneities of the planet to time- variations in disk-integrated 
brightness measured by a distant observer. In particular, we consider changes in disk-integrated thermal flux due to spatial 
inhomogeneities in thermal emission, and variations in disk-integrated reflectance due to spatial inhomogeneities in albedo. 

Russell (1906) showed how time-resolved photometry of reflected light constrains the albedo markings, shape, and ro- 
tational orientation of minor Solar System bodies observed near opposition (full phase). In addition to revealing brightness 
markings on stars and planets, rotational and orbital phase variations likewise have the potential to constrain rotational 
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and orbital inclination. Possible applications include: the thermal phase variations of non-transiting hot Jupiters might hint 
at their orbital inclination, breaking the Afsini degeneracy and allowing for improved mass estimates; the rotational phase 
variations of a transiting planet's host star may be sufficient to infer its rotational inclination (a.k.a. stellar obliquity), a useful 
discriminator between planet migration scenarios (Winn et al. 2005); the rotational photometric variations of a directly- imaged 
planet might encode information about its rotational inclination which — when combined with the astrometrically inferred 
orbital inclination — provides an estimate of planetary obliquity, telling us about planet formation (Tremaine 1991). 

1.2 Forward vs. Inverse Problem 

Inferring the properties of a star or planet based on its disk-integrated brightness is an inverse problem, as opposed to the 
forward problem of predicting the photometry of an object based on its properties. There are two broad approaches to solving 
this inverse problem: 1) develop a forward model with unknown parameters, then use a x'^-minimizer or Markov Chain Monte 
Carlo (MCMC) to solve for the parameters, or 2) approximate the forward problem by an analytic transformation that has 
an explicit inverse. The advantages of the former method are that the forward model may be arbitrarily sophisticated, and 
parameter uncertainties can be readily estimated; the advantages of the latter are speed and the ability to explicitly identify 
degeneracies. A typical compromise is to adopt a simple linear approximation to the forward model, but use a MCMC or 
similar routine to obtain best-fit parameters and uncertainties. 

In this work we seek an analytic solution to the forward problem that would lend itself to an analytic inverse. This entails 
determining the basis lightcurves for a set of complete and orthonormal basis maps. We approximate the forward problem as 
linear in the unknown planet map, M{6,cj>): 

F{t)= i K{e,(j>,t)M{6,(j))dQ., (1) 

where F{t) is the observed fiux, K{6, <j}, t) is the kernel, 6 and <j} are planetary co-latitude and longitude, respectively. The 
inverse problem, solving for M given K and F, is a Fredholm integral equation of the first kind and is non-trivial (Aster, 
Borchers & Thurber 2013). 

As we will see below, the kernel is non-negative and unimodal. It is therefore tempting to think of (1) as a convolution, 
and the inverse problem as a deconvolution. For thermal lightcurves, K has a fixed shape and this description is formally 
correct; in other cases it is merely a useful analogy. 

If the orientation of the planet or star is not known a priori, then the problem can be expressed as 

F{t) = <£ K{G, 9, (j), t)M{e, cj>)dQ., (2) 

where G represents the unknown geometry (inclination or obliquity). The object is then to solve for G and M{d, (j>) knowing F 
and the form, of K. It has been demonstrated in numerical experiments, for example, that one can simultaneously constrain a 
planet's two-dimensional albedo map, obliquity and obliquity phase (Kawahara & Fujii 2010, 2011; Fujii & Kawahara 2012). 
In the current work we seek an analytic solution to the forward problem, and hence a rapid and elegant approach to the 
inverse problem. Whenever possible, we solve the integrals analytically by hand and/or with Mathematica. When symbolic 
solutions are too messy to have intuitive value, we use IDL to compute and plot numerical integrals. In §2 we describe our 
chosen basis maps and lightcurves. We tackle thermal lightcurves in §3 before addressing the more complex case of reflected 
lightcurves in §4. In both of those sections we begin by describing our model assumptions for the forward problem, then 
present solutions to special cases before moving on to the general solution. Some of the mathematical details are relegated to 
the Appendix in order to keep the main body of the manuscript accessible. We discuss possible applications of this work in 
§5. 



2 HARMONIC LIGHTCURVES 

In order to develop an analytic solution to (1), it is necessary to express the planetary map analytically. In general this is 
done by decomposing M using an orthonormal basis. The obvious basis maps for a spherical planet are spherical harmonics. 
Any continuous, static albedo map, M{9,(j>), niay be decomposed as 

Mie, 0) = E E crrne, 0), (3) ^" = i / ^(^' ^)^™(^' ^)'^"- w 

1=0 m = -l 

The real spherical harmonics are given by 

Y"^(ff ^^ ^ / ^r^!m(cos 8) cos(m0) if m > 

' '^"'W \ ivl'"lp;|„|(cos6l)sin(im|0) if m < 0, ^ ' 

where Pim is the associated Legendre polynomial without the Condon-Shortley phase, (—1)™. 
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We adopt the geodesy normalization (unit power) for real spherical harmonics, 

[ V — {r+7nT< — It (> 0, 

The goal of this project is to determine the lightcurve signatures of spherical harmonics, or harmonic lightcurves, 

Fr{t)^{,K{eA,t)Yr{eA)dn. (8) 



It is perfectly equivalent to think of this as decomposing the kernel into spherical harmonics. 
Since (1) is linear in M, the lightcurve for an arbitrary albedo map is simply 

oc I 

^(*) = E E cTFrit). (9) 

The dream of exo-cartography is to take the inner product of an observed lightcurve with known harmonic lightcurves 
to determine map coefflcients: 

cr = ^ f Fr{t)F{t)dt, (10) 

where P is the period of the lightcurve. A map of the planet could then be reconstructed using (3). In this paper we will 
encounter two problems with this strategy. First of all, some of the harmonic lightcurves are flat lines [F = 0). In other 
words, the transformation has a nuUspace.^ Secondly, even non-zero harmonic lightcurves are sometimes proportional to each 
other. This is not surprising, since linear transformations need not preserve angles: planetary maps that are orthogonal are 
not necessarily transformed to lightcurves that are orthogonal. The bottom line is that exo-cartography suffers from formal 
degeneracies, even in the limit of perfect observations. 



3 THERMAL LIGHTCURVES 

3.1 Model Formalism 

We assume a static brightness map, diffuse thermal emission, and neglect limb-darkening. The requirement of a static map 
depends on context. For mapping star spots or patchy clouds on a brown dwarf, the rotation rate is the relevant timescale. 
When mapping the diurnal heating pattern of a planet, on the other hand, one requires stability on the orbital period (for 
more about the various sources of planetary thermal variability see Cowan et al. 2012c). 

The flux, F, in this case is the disk-integrated thermal flux from the planet. The kernel is proportional to the visibility of 
a given region of the planet at time t: K{9,4>,t) = i\/(^, 0, i), where the visibility, V, is unity at the sub-observer location, 
drops as the cosine of the angle from the sub-observer location, 70, and is zero on the far side of the planet: 

V{6, <j},t) = max[cos7o, 0] — max[sinSsin6oCos((^ — <^o) + cos ^ cos ^o, 0], (11) 

where 60 and (po are the sub-observer latitude and longitude, respectively. The piece-wise defined nature of the kernel leads 
to much of the difhculty in solving the forward problem analytically. 

The entire time-dependence of the forward problem comes in through the sub-observer position. In the absence of 
precession, the sub-observer latitude is constant, 9o{t) = 9o. The sub-observer longitude decreases linearly with time (we 
define longitude increasing to the East, with the planet rotating from West to East): <j}o{t) ~ (t>o{0) — tJroti, where iUiot is the 
rotational angular frequency in an inertial frame (e.g., uj^ot = 27r/23.93 hr~^ for Earth). 

Given a thermal map of the planet, the harmonic lightcurves are given by 

Frit) = ^fv{e,c^,t)Yr{e,<t>)dn. (12) 

Integrating the piece- wise defined kernel over the entire sphere is equivalent to integrating the non-zero part of the kernel, 
Knz{9,(j>,t) = -^ (sin 61 sin 6*0 cos (0 — <^o) + cosScos^o), over the visible hemisphere. The limits of integration are then defined 
by the limb, the locus of points with jo ~ 7r/2. From (11), the limb satisfies: 

taneiimb = 7 ^ 7- —^, (13) 

tanfyocos(<p — (po) 

as shown in Figure 1. 

For a planet viewed equator-on {6o = vr/2) the kernel simplifies to: Knz{d,(t>,t) ~ — sin 8 cos{(f) — (f>o)- 

Tiic term "kernel" is often used intcrciiangcably witii "nullspace" in matiiematical physics, but we escliew that terminology here 
because "kernel" already has a central role in convolutions. 
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Kernel 



Kernel 



a> 71/2 



371/4 




a> 7t/2 



371/4 




Figure 1. The visibility (red contours) of a planet seen by an equatorial (left) or northern observer (right). The gray-scale shows the 
kernel of the convolution, which in this case is simply the rescaled visibility. The sub-observer location is denoted by a red cross. The 
limb is shown in each case with green dotted lines. For either viewing geometry, the limb is a great circle and the non-zero portion of 
the kernel is a hemisphere. 

3.2 Equator-On Thermal Lightcurve 

We first consider a planet viewed equator-on (left panel of Figure 1), which allows us to separate (12) into two single integrals: 



iV," 



F^it) = i^ / y/T^^Pi^{x)d: 



0o + f 

■X I cos((j} — (f>o) cos^ni' 



(14) 



where we have made the change of coordinates x — cosO. Note that we have given the example for a cosine Y/" {m ^ 0), but 
the sine instance {m < 0) can be trivially obtained by replacing m — >■ \m\ and cos{m(j)) — > sm{\m\<j)). 

The time-dependence can be factored out of the integral by making the change of variables $ = — 0o, using a trigono- 
metric identity, and noting that one of the two resulting integrals is zero (or simply using complex exponentials): 



7?™(t) = _!_cos(m(?!)o) / Vl-x'2Pirnix)dx / cos($) cos(m$)d$, (15) 

where the product of integrals is now solely a function of I and m. Recalling that (j>o oc t, Fourier analysis implies J^ F{^F^dt = 
for /i 7^ 771, and justifies the use of sinusoidal basis maps by Cowan & Agol (2008). 

The integral of the Associate Legendre Polynomial is simplified by using a recurrence relation, then solved directly 
following Jepsen et al. (1955), as described in Appendix A: 



\/l - x^Pim(x)dx = 



{21 -f 1) 



•■i+i 



m+l] 






(16) 



The integral of an associated Legendre polynomial over the interval x £ [—1,1], represented here as _R™, depends on the parity 
of the spherical harmonic. If / and m are even, then I — 1, ? + 1, and in + 1 are odd; if I and m axe odd, then 1 — 1, I + 1, and 
m -I- 1 are even; if / -)- m is odd, then so are {I - 1) + {m + 1) and {I + I) + {m + 1). Lastly, i?™^^ = ^I+t^ fo'^ odd I > 1, 
producing a nuUspace, F™ — 0. 
The "l>-integral yields: 



cos $ cos(m<l>)d$ = 



if Iml 



1 



■cos(^) if|ml/l, 



(17) 



which recovers the Cowan & Agol (2008) result of zero phase signature for odd \m\ > 1, since cos (^) is zero in those cases. 
Physically, the brightness inhomogeneities cancel each other in the disk-integrated case. Mathematically, Y{^ with odd |77i| > 1 
are in the nuUspace of the convolution. 

Combining (16) and (17) yields the full solution, 



r 1 



Frit) 



V3 ""= '^' 
2{_l)"»/2 



if ; = o 

a I = 1 and m — 1 



2(l~m)\ 



^ii-m'^) \l (2i+i)(i+m)! LKlt^ (o^d) - i?™+^(odd)] cos(m(^o) if I and m are even 



(18) 







otherwise, 



where _R,'"(odd) is given in Appendix A. The nuUspace is the union of odd m > 1 ($- integral goes to zero) and odd I > 1 
(a;- integral goes to zero). 
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The first few non-zero harmonic lightcurves are: 



FSit) 


= 1 


Fi\t) 


2 A 

= —j=COS<po 

v3 


F^it) 


^/T0 

8 


Fi{t) 


— —-— cos{2(po 
8 



-Ht)- 


3V2 
64 


(23) 


"iit) - 


- 32 cos(20o) 


(24) 


4 (i) -- 


64 '=°^('*''^°)' 


(25) 



(19) 
(20) 

(21) 

(22) 

where again we stress that the sine harmonic lightcurves {m < 0) can be trivially obtained by replacing m -^ \m\ and 
cos — > sin. 

3.3 Inclined Thermal Lightcurve 

For a general (non-equatorial) observer, one of the poles is visible and the other is invisible, rather than both being on the 
limb (right panel of Figure 1). The (/)-integral therefore runs from to 27r, while the S- integral runs from to Siimb(0) if 
the north pole is visible, or from ^iimb(0) to tv if the south pole is visible. With no loss of generality we consider a northern 
observer: 



Frit) =^ 

TV 



2ir 

sin6^o 



pZTV pi 

/ cos((jf) — (jio) cos(m(j!)) / \/l — x^Pim{x)dx 



/'27r nl 

cosOo I cos{m(f>) / xPim{x)dxd(p 

JO "'^limb 



(26) 



where ajumb = cos6liimb- 

As with the equatorial geometry, the time-dependence may be factored out of the integral using complex exponentials: 

Arm r ^27r ^1 



sm Ho 



f- ZTT pL 

/ cos$cos(m$) / \/l — x^Pim{x)dxd^ 
■^0 J^nmb 



TV 

(27) 



-f cos^o / cos(m<I>) / xPim(x)dxd<^ 

■^0 J X\im\, 



where (13) dictates that 



— tan So cos <& ,„„. 

a^iimb = , „ „ ■ (28) 

Vl + tan2 Oo cos2 $ 

Since ajnmb is a function of "1>, we cannot separate the x and $ integrals as in the equatorial case. We instead denote the 
two double integrals as 

P'2tt pi 

Ar{Oo)= cos<&cos(m$) / Vl - x'^Pim{x)dxd$, (29) 

p2n pi 

BTido) = / cos(m$) / xPi^{x)dxd^. (30) 

Numerically, we find that A = B = Q li and only if I is odd and greater than 1, therefore all thermal harmonic light curves 
are zero for odd I > 1 (we demonstrate this analytically in Appendix B). The nuUspace for an inclined geometry is therefore 
more limited than the equator-on case: there are non-zero harmonic lightcurves with odd m> 1 (provided that I is even, e.g. 
Yi shown in Figure 2). 

As with the equatorial case, one can use recurrence relations to express the a;- integrals in terms of sums of simple integrals 
of Pr{x). But while there are recurrence relations for the indefinite integral j Pr{x)dx (DiDonato 1982), we were unable to 
develop such a relation for the harmonic lightcurves, Fl^it)- Instead, we use the brute force approach of substituting specific 
Pi"{x) into (27) and analytically solving the double integral. 

We first solve the indefinite integral over x, which we evaluate at the limits of integration. The resulting integrands for 
the "l>-integrals include trigonometric functions with singularities. Since the sign of these functions can change on either side 
of the singularities, the integral over $ exhibits unphysical jumps. Fortunately, the jumps occur at predictable fractions of tv, 
and their amplitude is a tractable function of 9o, so the definite integral can be manually corrected. The resulting low-order 
non-zero harmonic lightcurves are: 

^ Although these jumps are undesirable, the derivative of the resulting curve is smooth, so Mathematica is indeed returning a valid 
anti-derivative. 
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Figure 2. The Y^'* brightness map is in tlic nullspacc of tlie edge-on thermal lightcurve, but not of the inchned lightcurve. It is the 
lowest-order map that can produce power at odd m > 1 in a thermal lightcurve. 
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The sine harmonic hghtcurves (m < 0) can be triviaUy obtained by the substitution cos(m 



(37) 
(38) 

(39) 
(40) 

(41) 



-> sm{\m\<f>o)- 

In Figure 3 we show the amplitude of low-order harmonic hghtcurves as a function of sub-observer latitude, 60 {Oo = 
for pole-on or face-on rotation; 60 ~ 7r/2 for equator-on or edge-on rotation). The nullspacc of the convolution, and more 
generally its eigenvalues, is a function of 60 ■ Obviously a pole-on object exhibits no rotational variability (left side of the plot). 
Moreover, an object exhibiting power in m = 3 is neither pole-on nor equator-on: the amplitude of F4 exhibits a clear peak 
at 60 ~ tt/3. 

Formally, there are an infinite rmmber of harmonic maps that contribute to the lightcurve power at a given m. However, 
the low-pass nature of the convolution ensures that most of the power at a given m comes from the lowest-order F,™. If the 
intrinsic power at Y/ was known a priori and there was no other map contributing m = 3 power, then one might hope to 
estimate rotational inclination simply by measuring the m = 3 power present in the lightcurve. But there would still be a 
two-way degeneracy because 60 is double- values for a given amplitude (green line in Figure 3). 



4 REFLECTED LIGHTCURVES 

4.1 Model Formalism 

We assume a spherical planet with a static albedo map on a circular orbit. If one is performing rotational mapping (Cowan 
et al. 2009) then the map only needs to be constant over a single rotation; or equivalently the recovered map is a diurnal 
average. For rotational-|-orbital mapping ("spin-orbit tomography;" Fujii & Kawahara 2012), the map is assumed to be static 
over an entire planetary orbit (e.g., Hasinoff et al. 2011). 

The reflection is treated as diffuse (Lambertian). Note, however, that real surfaces can exhibit specular reflection (Williams 
& Gaidos 2008) and atmospheres can exhibit Rayleigh or Mie scattering (Robinson et al. 2010). The albedo map should be 
thought of as a top-of-atmosphere planetary albedo. 

Albedo, and hence all of the observables, will in general be a function of wavelength, but since we consider only scattered 
light, there is no mixing of wavelengths. Our results can be generalized to any number of wavebands (Kawahara & Fujii 2011; 
Fujii & Kawahara 2012) or arbitrary linear combinations of wavebands (Cowan et al. 2009, 2011; Kawahara & Fujii 2010; 
Cowan & Strait 2013). 

For a uniform planet the resulting phase variations under the assumption of diffuse reflection is the well-known Lambert 
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Figure 3. The relative amplitude of harmonic thermal lightcurves as a function of sub-observer latitude (or, equivalently, rotational 
inclination). Line style denotes I, color denotes m. The nuUspace of the convolution is a function of do, allowing for qualitative estimates 
of rotational inclination based on the Fourier spectrum of observed lightcurves. 



phase (Russell 1916). In this section we derive higher moments of the lightcurve. This exercise is complementary to computing 
the phase variations of uniform but non-Lambertian planets (e.g., Tousey 1957; Madhusudhan & Burrows 2012). Real planets 
are neither uniform nor Lambertian. 

The reflected-light kernel is K{9,(f>,t) = -V{6,cj>,t)I{6,(j),t), where the visibility, V is defined as in (11), while the 
illumination, I, is unity at the sub-stellar location, drops as the cosine of the angle from the sub-stellar location, 7s, and is 
zero on the night-side of the planet: 



I{6,(j),t) = max[cos7s,0] — maxfsinSsin^a cos( 



-I- cos cos ( 



,0], 



(42) 



and 9s and (ps are the sub-stellar latitude and longitude, respectively. Note that the kernel is proportional to the normalized 
weight, W , but is not divided by the Lambert phase function (Cowan et al. 2011). 
The reflected-light forward problem is therefore: 



^"(i) = ^<f v(e, <j,, t)i{e, <^, t)Yr(e, ^)dn, 



(43) 



where F in this case is the reflectance of the planet, i.e. the planet/star contrast ratio after accounting for the radius and 
semi-major axis of the planet, which are sometimes not of immediate interest. 

Noting that the orbital phase, a € [0, tt] (a = at full phase; a = tt at new phase), is simply the angular distance between 
the sub-observer and sub-stellar points, it may be expressed as 



cos a = sm do sm Os cos( 



-f cos 6o cos 9 s — cos ^ sin i, 



(44) 



where the constant i G [0, 7r/2] is the orbital inclination with respect to the celestial plane {i — for a face-on orbit; i = 7r/2 
for edge-on), and ^ — ^(0) + iUoibt is the planet's orbital location (^ G [0, 27r]; ^ = at superior conjunction; ^ = tt at inferior 
conjunction). 

As with the thermal emission problem, the crux stems from the piece-wise defined kernel, or equivalently, the limits of 
integration for the non-zero portion of the kernel. The reflected lightcurve calculation is harder than the thermal lightcurve 
because the non-zero region of the integral is a lune rather than a hemisphere. 

The analytic expressions for the limits of integration allow us to drop the piecewise-defined version of the visibility and 
illumination functions. If one only considers those regions where both illumination and visibility are greater than zero (ie: 
those within the limits of integration), then the kernel is: 



TrKnz{9,(l>,t) = (sin S sin So cos ( 



cos 6 cos Oo ) (sin 6 sin 6s cos ( 



' cos y cos fs ) 



(45) 



In the current study, we only consider the tidally locked configuration, which is likely to be relevant for the current 
crop of hot Jupiters as well as temperate planets orbiting low-mass stars. In this case, the sub-stellar location is equatorial, 
ds = 7r/2, 9o — i, and one can place the prime meridian at the sub-stellar meridian with no loss of generality {(jjs — 0): 
T^Knz {6,(f>,t) — sin i sin^ 9 cos (f> cos{(f> — (f>o) + cos i sin 9 cos 9 cos (f>. 

If the planet also orbits edge-on, (fis ~ 0, then: ■KKnz{9, 4>, t) = sin^ 9 cos cos (0 — (f>o)- 
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Figure 4. The visibility (red contours) and illumination (yellow contours) of a tidally-Iocked planet seen by an equatorial (left) and 
northern (right) observer. The limb and terminator are shown in each case with green dotted lines. The sub-observer location is denoted 
by a red cross, the sub-stellar location by a yellow astrerisk. The limb and terminator are both great circles, making the non-zero portion 
of the kernel a lune. The maximum value of the kernel, and its integral, are a function of the angular distance between the sub-observer 
and sub-stellar locations, i.e.: the planet's orbital phase. 



4.2 Tidally-Locked, Edge-On Reflected Lightcurve 

The combination of tidally-Iocked rotation and edge-on orbit dictates that the 6'-integral runs from to tt, so the double 
integral can be split into two single integrals: 



Fr[t) = 



7=" 7-1 



2 /"^^ 

(1 — X )Pim{x)dx / COS 4> cos{(f> ~ (f>o) cos{m(j>)d(j>, 



(46) 



where (j>i — max[— 7r/2, </>(, — 7r/2], and (j)2 = min[7r/2, </>(, -|- n/2] (left panel of Figure 4), and we have again shown the 
cosine case. The x-integral is merely a scalar, so -F™ oc Fi^, as for the edge-on thermal case. The time-dependence cannot be 
factored out of the (^-integral, however, so the negative and positive m lightcurves are not trivially related by the substitution 
cos(7n0) — ^ sin(|7Ti|(^) (cf. the m = —1 and m — 1 cases of Equation 47). 

As with the edge-on thermal lightcurve, the integral over x may be solved for arbitrary I and m by using recurrence 
relations for Associated Legendre Polynomials and the Jepsen et al. (1955) solution to their definite integral. In the present 
case, this leads to terribly cumbersome expressions, so we instead adopt the brute force approach of analytically integrating 
the x-integral for specific I and m. 

The (^-integral must be solved separately in two cases: (po < and (f)o 1^ for the first and second halves of the planet's 
orbit, respectively. The two cases may be stitched together by noting that a = \4>o\'- 



,'1 



■d4>^ < 



^|i^il^((m + 2)sin(a- 

i sin 00 (tt — a + sin a cos a) 

I sin 00(1 + cos a) 

i ( sin Q + (tt — a) cos a J 

|cos^(<^„/2) 

J cos <^o (tt — a + sin ct cos a) 

^fi=|^((m + 2)sin(a- 



rag 
2 






m-K 
2 . 



+ (m — 2) sin (a + 



ma 
2 



mn 
2 



)) 



2 . 



+ {m — 2) sin (a + 



ma 
2 



mn \ 
2 ) 



if m < — 2 
if m = ~2 

if 771 = —1 

if m, = 

if 771 = 1 
if 771 = 2 

if 771 > 2. 



(47) 



The first few non-zero harmonic lightcurves are: 
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Figure 5. Harmonic reflected lightcurves for a tidally-locked planet on an edge-on orbit. Line style denotes I, color denotes m. 
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(50) 
(51) 
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32 V 6 
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^3(i)=^\/^cos*(0o/2) 



Fi W = Y^ Y Y (4 cos 0O + 5 cos(20o) - cos(40o) 
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'I 
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sm a cos a 



TT V35 sm ( 



sin acos(20o)- 



(54) 
(55) 
(56) 
(57) 
(58) 
(59) 



7rV35 

A uniform map, M(Q, 0) = 1, produces the Lambert phase function, as expected. The nuUspace of the convolution includes 
odd I + m (equivalently, odd I — m; i.e., an odd number of nodes in the meridional direction), for which the contribution 
from the northern and southern hemispheres cancel perfectly. Unlike the edge-on thermal case, and counter-intuitively, purely 
meridional maps, such as Y2 {6, 4>), are not in the nuUspace. 

The harmonic lightcurves are shown in Figure 5. As expected, the negative and positive m are not related by a simple 
phase shift. Significantly, harmonic lightcurves with the same m are proportional to each other. Therefore, although the 
nuUspace of this convolution is more limited than for thermal lightcurves, there are still severe formal degeneracies in the 
inverse problem. 



4.3 Tidally-Locked, Inclined Reflected Lightcurve 

In the more general case of a tidally-locked planet on an inclined orbit, the expression for harmonic lightcurves becomes: 



Fr{t) = -^i sini 



r/2 



COS COS 



w/2 



{1- X )Plm{x)dx 



d<t) 



+ cosi 



r/2 



7r/2 



im4> 



I x\Jl - x'^Pi^{x)dx 



(60) 



where a;iimb is given by (28) and the meridional limits of integration implicitly assume a northern observer (right panel of 
Figure 4). For a southern observer the limits would be [—1, arumb]- The limits of integration for the 0-integral no longer span 
27r radians and therefore do not lend themselves to the Fourier strategy used for the inclined thermal lightcurves."^ Instead, 
we solve the integrals numerically and plot the solutions in Figure 6. 

^ Moreover, the locations of the jumps in the integrand of the 0-integral are not as predictable as they are in the inclined thermal 
lightcurve case. 
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Figure 6. Harmonic reflected lightcurves for a tidally locked planet on an inclined orbit (i = So = ^)- Line style denotes I. color denotes 



Prom numerical integration, we determine that there is no nullspace for the inclined reflected lightcurves up to / = 4, nor 
do the lightcurves depend simply on do- In other words, one cannot simply scale the edge-on solutions of Figure 5 by sin i. 
For example, the shape of F^ remains fixed but with the addition of a constant offset. This means that the Y^ component 
of the map has the effect of decreasing the reflectance of the planet at full phase and increasing it at crescent phases. This 
is essentially a restatement of the latitude-albedo effect noted by Cowan et al. (20f 2b): reflective poles make a low-obliquity 
planet appear abnormally bright at crescent phases. 

5 DISCUSSION 

5.1 Thermal Phase Variations of Short Period Planets 

It is significant that there are different nuUspaces for the edge-on and inclined thermal lightcurves. One can't simply scale 
the edge-on solution of Cowan & Agol (2008) by sini. For transiting systems, the orbital inclination is nearly 90° and is well 
measured by the transit morphology (e.g., Charbonneau et al. 2000). This means that while the amplitude of odd modes in 
an observed thermal lightcurve is likely an order of magnitude weaker than for adjacent even modes, a measurement of these 
modes could be readily converted into an estimate of the N-S asymmetry of the planet. 

Moreover, the amplitude of odd modes in the observed lightcurve of a non-transiting planet provides a qualitative 
estimate of orbital inclination, since the planet may be assumed to have zero-obliquity. The original motivation for measuring 
thermal phase variations of hot Jupiters was to estimate orbital inclination in order to break the M sin i degeneracy (Agol 

6 Charbonneau 2005). That proposal was based on the presumption that all hot Jupiters would have the same day-night 
temperature contrast (m = I amplitude). While that assumption is demonstrably wrong (Cowan et al. 2007), the inclination- 
dependance of thermal phase variation Fourier spectra offers an opportunity to constrain orbital inclination after all. 

The best-characterized hot Jupiter, HD 189733b, only has empirical constraints up to m = Z = 2 at 3.6 and 4.5 micron 
(Knutson et al. 2012). The amplitude of the m = 2 component of the phase variations is 5% of that at m = 1 for both 
wavebands, indicating that the m = 2 map has only 12% the amplitude of the m — 1 map (following Cowan & Agol 2008). 
This is not surprising since diurnal forcing primarily excites the m = 1 mode. The 8 micron 2D map of HD 189733b only 
constrains lightcurves to Z = 1 (Majeau et al. 2012), so the strength ofm<l modes, in particular the N-S asymetry responsible 
for odd lightcurve modes, is currently unknown empirically (and debated theoretically: Cho et al. 2003; Cooper & Showman 
2005). Aside from the intrinsically weak signal, possible complications include limb-darkening, which slightly modifies the 
convolution kernel (Cowan & Agol 2008), the presence of eccentricity seasons (Lewis et al. 2013), or the contamination of 
m = 2 modes by ellipsoidal variations (Cowan et al. 2012a). 

5.2 Rotational Modulation of Stars and Brown Dwarfs 

Likewise, the Fourier spectra of spotty stars (Lanza et al. 2009) and brown dwarfs (Artigau et al. 2009; Radigan et al. 2012) 
might hint at their rotational inclination. This is unsurprising given that star spot modeling can yield inclination estimates 
good to tens of degrees (Walker et al. 2007). However, aside from being faster than traditional spot modeling, our Fourier 
approach suggests why star spot modeling is able to constrain inclination despite the myriad degeneracies (Dorren 1987; 
Kipping 2012). It should be feasible, by the same token, to measure the presence of odd modes in Spitzer Space Telescope 
lightcurves of spotty brown dwarfs to obtain a qualitative measurement of their inclination. 
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Figure 7. Fourier spectrum for a star of unit intensity witii a square spot of size tt/IO, zero intensity, and centered on 6 = vr/S, = vr/S. 
Each color corresponds to the power spectrum as measured by an observer at a different inclination. The spectra have been normalized 
to their m = 1 amplitude. An edge-on rotation produces no power in the odd modes (the non-zero values on the plot are due to numerical 
errors). The plot looks qualitatively similar for different spot sizes and locations, provided the spot is not equatorial, which kills m j^ I 
modes in the map, and hence odd modes in the lightcurve. 



The stellar/brown dwarf inclination inverse problem is more favorable than extracting the orbital inclination of hot 
Jupiters because the signal-to-noise ratio for the variations are much greater (~ 1% rather than ~ 0.1%), and the intrinsic 
power spectrum of the thermal map does not drop precipitously with m. The limiting case of a single 5-function spot has 
a flat Fourier spectrum. This makes Fourier analysis an inefficient means of mapping star spots, but potentially useful for 
constraining rotational inclination. Lumme et al. (1990) used the power spectra of reflected lightcurves to extract the rotational 
orientation of asteroids, independent of shape and albedo markings. Our work extends this approach to thermal lightcurves. 

In Figure 7 we show the power spectra for a rotating star at different inclinations. For stars that are rotating nearly 
edge-on {9o = 90-80° , denoted by solid lines) the power at odd m is a very strong function of orientation. If one reasonably 
assumes that a star's intrinsic intensity map should not be deficient in odd modes, then a factor-of-two measurement of odd 
amplitudes constrains the rotational inclination to within a degree. While the a priori odds of a star being within 10° of 
edge-on are not great, they increase substantially for stars known to host transiting planets (Sanchis-Ojeda et al. 2012). At 
larger inclinations (denoted by dotted lines), the orientation of the star still leaves an imprint on it's Fourier spectrum, but 
one would need a strong prior on the intrinsic power spectrum of the stellar intensity map. It is not clear whether a Fourier 
approach would be advantageous to the usual star spot modeling in such a case. 

Fourier analysis of rotating inhomogeneous bodies is analogous — conceptually and mathematically — to estimating stellar 
inclinations with asteroseismology (Gizon & Solanki 2003). Although our proposed method requires high-precision photometry, 
the observational cadence could be considerably lower than for asteroseismology, so many Kepler target stars may be amenable 
to such an analysis, without the need for v sin i measurements from high-resolution stellar spectra (e.g., Hirano et al. 2012). As 
with any star-spot based inclination estimate, problems include differential rotation of star spots, or their formation/dissipation 
on rotational timescales. If the brightness markings evolve on timescales longer than the rotational period, it should still be 
possible to Fourier decompose each rotation separately to perform the spectral power analysis. Limb-darkening is again a 
source of systematic uncertainty/error. The intrinsic power in the Y^ map will depend on the latitude of star-spots, with 
equatorial spots providing no power at this harmonic. But star spots are, if anything, less likely to be at the equator (Baumann 
et al. 2004). It is not known whether the clear/cloudy regions on L/T transition dwarfs have preferred latitudes. 

5.3 Obliquity of Directly-Imaged Planets 

Directly imaged jovian planets emit thermal radiation because they are still young. Although no rotational variability has 
yet been reported for these objects, they are probably cloudy (e.g., Madhusudhan et al. 2011) and may exhibit the same 
sort of variability as brown dwarfs. Merely detecting this rotational modulation would be a technical feat with important 
implications for giant planet formation. Further down the road, however, we might hope to measure their rotational lightcurves 
with sufficient precision to constrain their rotational inclination with respect to our line of sight. Unlike hot Jupiters, directly 
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imaged Jovian planets are not tidally locked, but their orbital inclination may be estimated observationally. The Fourier 
spectrum of such a rotational lightcurve would therefore put a joint constraint on the planet's obliquity and equinox via: 

cos^o = sin O cos ^e sini + cos O cos i, (61) 

where is the planetary obliquity and ^e is the angular location of northern summer solstice with respect to superior 
conjunction. 

This is not the only method that has been proposed to measure the obliquity of directly-imaged planets. Full-orbit 
thermal phase curves of mature planets, for which insolation rather than internal heat dominates the power budget, might 
in principle betray the obliquity of a planet, but the inverse problem is complicated by orbital eccentricity, diurnal heating, 
and the details of heat storage and transport (Gaidos & Williams 2004; Cowan et al. 2012c). Reflected full-orbit lightcurves 
of directly-imaged planets have also been demonstrated to convey information about planetary obliquity (Fujii & Kawahara 
2012). These methods have the advantage that they can break the degeneracy between obliquity and its orientation, but they 
require observations spanning the planetary orbit rather than its rotation, and therefore may not be practical for long-period 
planets. 

5.4 Reflected Phase Variations of Planets 

The fact that F2 oc Fq for edge-on reflected phase curves has important implications for the retreival of albedo from planetary 
phase variations. The Y2 map corresponds to a zonally uniform planet with bright poles and has the effect of reducing the 
amplitude of zeroth-order (Lambertian) phase variations. Consider the worst-case of a diffusely-reflecting planet with albedo 
map M — \YS ^ J— Y2, which has an albedo of unity at both poles, zero at the equator, and a mean albedo of (M) = i (this 
is qualitatively similar to the albedo map of a low-obliquity planet with polar snow/ice). The Bond albedo of the planet is 
A = ^A'i + j^A9. = i (Appendix C). The reflected lightcurve for the planet, however, is F{t) = iF,?(t)-|-— 2=F^'(t) = ^F^\t). 
In other words, an observer would see a planet exhibiting perfectly Lambertian phase variations and would infer a planetary 
albedo of 20%. The sub-observer and sub-stellar latitudes are both equatorial, but the albedo estimate differs from the actual 
Bond albedo by 0.05/0.20 — 25% (and differs from the mean albedo by 67%). One should therefore be wary of estimating 
an exoplanet's Bond albedo, even if the planet orbits edge-on, is tidally locked, and exhibits Lambertian phase variations. 
Precision is no guarantee of accuracy. 

The inclined tidally-locked reflected lightcurves exhibit no nuUspace up to / = 4, but the non-orthogonality of harmonic 
lightcurves leads to formal degeneracies for mapping planets, even in the limit of noiseless data. It remains to be seen to 
what extent such degeneracies affect the spin-orbit exo-cartography of Fujii & Kawahara (2012). Although we have only 
derived reflected lightcurves for tidally-locked planets, our results regarding meridional albedo markings affecting reffected 
phase variations should apply for non-locked planets, provided one averages over the rotational variation, as well as non-zero 
obliquity (e.g., the numerical simulations of Cowan et al. 2012b). 

By deriving the reflected lightcurves exhibited by spherical harmonic maps, Russell (1906) demonstrated the existence 
of a nuUspace for albedo markings on a distant body viewed at opposition: the odd harmonics. Researchers have dealt with 
this formal degeneracy by imposing positivity of albedo at every location on the body, adopting a limited number of discrete 
spots (Marcialis 1988), and/or using only two albedo values (light and dark regions; Lacis & Fix 1972). In this paper we have 
shown that such degeneracies are generically present in exo-cartographic inverse problems. We have not attempted to resolve 
these degeneracies, since the specific strategies involved will likely depend on the application. The silver lining is that the 
inclination-dependence of the convolution may provide leverage for constraining the viewing geometry of unresolved systems. 
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APPENDIX A: INTEGRALS OF ASSOCIATED LEGENDRE POLYNOMIALS 

A technical crux of the analytic forward problem is solving definite integrals of Associated Legendre Polynomials, Pim- 

Al Recurrence Relation 

Recurrence relations allow us to relate functions of x including Associated Legendre Polynomials to combinations of simple 
Associated Legendre Polynomials. For example, we used 

Vl - x^Pi^ = ^y^ [Pi-1,^+1 - Pl+l,m+l] ■ (Al) 

A2 Definite Integrals on a; G [—1, 1] 

For the special cases where the limits of integration are x £ [—1,1] (edge-on, zero- obliquity reflected light; equator-on thermal 
emission), compact solutions (i.e., not involving sums) have been worked out by Jepsen et al. (1955). Those authors solved 
the deflnite integral of Pim{x) without the Condon-Shortley phase, precisely what we need in the current paper: 

[ ^r(even) = ,p3iy^§|0t=L if / and m are even 

Rr^l_^PUx)dx^l fir(odd) = - ,,.,^.,i,^-,';,+7>;!':;;-,i,+^),,i, if /and mare odd (A2) 

I if / -I- m is odd, 
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Physically, the northern and southern hemispheres have perfectly cancehng hghtcurves in the third case. 

APPENDIX B: WHY THERMAL HARMONIC LIGHT CURVES ARE ZERO FOR ODD I > 1 

The heart of this claim is the recognition that the calculation of (29) and (30) is formally equivalent to identifying the Fourier 
coefficients of the a;-integral. These Fourier coefficients will vanish unless the integral of the Legendre polynomial contains a 
term that has a periodicity that matches the term being pulled out by the Fourier integration, i.e. a term like cos m$ with 
the correct periodicity. We argue that this never happens for odd I > 1 {in fact, ^ = for I = 1 also). 
First consider our expression for ajumb, 

— tanSoCos<J> ,^^. 

a^iimb = , ^ „ . (Bl) 

Vl+tan2 6»oCos2$ 

As long as tan^o < 1 (i.e. 6o < vr/4) we can Taylor expand the denominator to obtain a power series in cos"!?, 

a^iimb ~ — tan Oo cos $ (1 — - tan^ do cos^ "^ + q ^an"' 9o cos** $ H J . (B2) 

Note that the power series within parentheses only contains even powers of cos $, and so when you multiply by the sole 
external term you obtain a power series containing only odd powers of cos $. Let's consider what this means for powers of 
a;iinib and one other function of it. 

The series associated to xf^^^, can be found by multiplying the two series for a^umb and because each term is an odd power 
of cos $ the product of two such terms will be an even power of cos <J>; thus the series for Xi^^^^ and all other even powers of 
a^iimb is completely made up of even powers of cos $. Similarly the series for x^^i^ and all odd powers of a^umb are completely 
made up of odd powers of cos<l?. Finally, note that 



f - xf,^^ = I = = f 1 - i tan^^ Oo cos^ $ + | tar/ 9, cos^ $ + ■■■), (B3) 

V 1 + tan^ 9c cos2 ^ \ 2 8 / 

which is a completely even series in cos$. Any power of this function will also be an even series in cos"l>. 

Any odd power of cos<l? can be expanded to obtain a series in terms of cosm<l? where m is necessarily odd, and similarly 
for even powers except that m is necessarily even and a constant term can also appear. 

To complete the argument the specific observations above have to be combined with several general observations about 
the associated Legendre polynomials. Below we will frequently encounter expressions containing (f — a;^)'"'^, which after 
Taylor expansion give an even power series in x. For this reason we will call these terms even, although it only applies after 
the expansion has been performed. From the Rodrigues' formula, 

(1 _ 2\m/2 im + l . , 

it follows that if Z + m is even then the Legendre polynomial is an even power series in x. On the other hand, if Z + m is odd 
then it is an odd power series in x. These basic observations about the associated Legendre polynomials are also valid for 
their indefinite integrals, the proof of which follows from the results of DiDonato (1982). Applying DiDonato's results we find 
that for the integral 

I sjl - x^Pi^{x)dx, (B5) 

if Z + 7n is even the result of the integration is an odd power series in x and vice versa. While for the integral, 

xPim{x)dx, (B6) 



if Z + m is even the result of the integration is an even power series in x and if it is odd then the result is an odd power series 
in X. (Note that for even I these integrals also give rise to arcsine functions and these must also be Taylor expanded for the 
statements to be valid.) 

If we adopt the complex exponential notation, the integrals of interest for the thermal light curves are 

Ar{eo)= f " cos <^e'""" f Vl-x^Pimix)dxd^ (B7) 

and 

Br{9o)= f "e™* / xPi^{x)dxd^. (B8) 

Putting aside the upper limits of the internal integrals for a moment and taking only the lower limits we can argue as follows: 
Focus on the A integrals first and note that for m ^ 0, 

cos<E>e™* = i('cos{(m-l)$} + cos{(?n+l)<I>}y (B9) 
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Thus viewing this as a Fourier decomposition the A integrals are the sum of the (m — 1) and {m + 1) coefficients. For odd I 
and even m, I + m is odd and so the Legendre polynomial integral in A gives rise to an even power series in x, which in turn 
gives rise to a sum over terms of the form cos (nx) where n is either zero or even. However, in this case (m ~ 1) and (m + 1) 
are odd and hence the integral is zero (there are no odd terms to pick out of the series). Similarly, for odd I and odd m, l + rn 
is even and we get a power series with only odd powers of x, which in turn only give rise to odd terms in cos (nx). And in 
this case (m — 1) and {m + 1) are even and once again the integral is zero. Thus, A™ = for odd I > 1. 

Chasing the even-odd parity oil, m and l + m similarly for the B integrals gives you the same result for I > 1, i.e. B™ = 
for odd I > 1. 

It is interesting to briefly comment on why these results don't apply to the case 1 = 1, and we also explain why the upper 
limits of integration don't contribute to the integrals considered above. For I = 1 and m = the argument above does apply 
to the A integral while the B integral is no longer formally equivalent to a Fourier series, there is no cosine term in the outer 
integral, and so F{' receives a contribution from the B integral. In the I = l,m = 1 case, the argument above applies to the 
B integral but not to the A integral which after using Eq. (B9) has a contribution with no cosine term. In both these special 
cases the upper limit of integration contributes because it is not integrated against a cosine term which when it is present 
cause the integration of the constant term to die. 

One final comment: the argument above started by assuming that 6o was in the range [0, 7r/4] in order that the Taylor 
expansion of the square-root made sense. However, the expansion can be continued to the range [7r/4, 7r/2] and the same 
argument run again, although the series is no longer convergent and hence is rather formal. In this manner, the conclusions 
about the null space for the northern observer can be carried over to the full range 9o € [0, 7r/2]. 

APPENDIX C: BOND ALBEDO OF A TIDALLY-LOCKED PLANET 

The Bond-albedo for a tidally-lockcd planet is time-invariable: 

A = -i A(e, 4>)i{e, 4>, t)dn. (ci) 

The contribution to the Bond albedo from a cosine harmonic map is 

AT = ^— \/l-x^Pi^(x)dx / cos(l>cos(m(j))d<j), (C2) 

^ i-1 i-f 

which is identical to (15) but with (po = 0. By analogy, the solution is simply 

( 1 if / = 



A7 



■y= if Z = 1 and \m\ = 1 



^(-l)"'''" . / 2(l-my. rpm + l/ , ,N _ „m+l 



^(1-m-^) V (2i+i)(i+m)! l^T+X (odd) - i?™^ (odd)] if I aud m are even 
otherwise. 



(C3) 



where i?["(odd) is given in Appendix B. The nuUspace is the union of odd \m\ > 1 ((^-integral goes to zero) and odd I > 1 
(a::-integral goes to zero). 



© 0000 RAS, MNRAS 000, 000-000 



